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Modeling interaction of relativistic and nonrelativistic 
winds in binary system PSR 1259-63/SS2883. 
I. Hydrodynamical limit 
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ABSTRACT 

In this paper, we present a detailed hydrodynamical study of the properties of the flow 
produced by the collision of a pulsar wind with the surrounding in a binary system. 
This work is the first attempt to simulate interaction of the ultrarelativistic flow (pul- 
sar wind) with the nonrelativistic stellar wind. Obtained results show that the wind 
collision could result in the formation of an " unclosed" (at spatial scales comparable 
to the binary system size) pulsar wind termination shock even when the stellar wind 
ram pressure exceeds significantly the pulsar wind kinetical pressure. Moreover, the 
post-shock flow propagates in a rather narrow region, with very high bulk Lorentz 
factor (7 ~ 100). This flow acceleration is related to adiabatical losses, which are 
purely hydrodynamical effects. Interestingly, in this particular case, no magnetic field 
is required for formation of the ultrarelativistic bulk outflow. The obtained results 
provide a new interpretation for the orbital variability of radio, X-ray and gamma-ray 
signals detected from binary pulsar system PSR 1259-63/SS2883. 

Key words: HD - shock waves - pulsars: binaries. 
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1 INTRODUCTION 

Pulsars lose their rotation energy through relativistic winds, 
the collision of which with the Interstellar Medium results 
in the formation of the Pulsar Wind Nebulae (regions of 
nonthermal synchrotron radiation of ultrarelativistic elec- 
trons accelerated at the termination of pulsar winds ( |Rees| 
fc Gunn|1974| |Kennel fc Coroniti|1994| l). The Crab Nebula 



is the most famous example of such an object. The recent 
X-ray and TeV gamma-ray observations ( jGaensler fc Slane| 
2006) show that this is a common phenomenon. 



A very interesting situation arises when a pulsar is lo- 
cated in a binary system. In this case the pulsar wind in- 
teracts with the wind from the companion star. This case, 
in particular, is realized in the binary system PSR1259- 
63/SS2883 which consists of a ~ 48ms pulsar in an elliptic 
orbit around a massive B2e optical companion ( |Johnsto"n| 
19941. The density and velocity of the stellar wind depend 



on the separation distance between two stars. Thus, the pro- 
cesses related to the interaction of two winds, in particu- 
lar particle acceleration and radiation proceed under essen- 
tially different physical conditions depending on the orbital 
phase. This makes this object a unique laboratory for the 
study of nonthermal processes in "on-line" regime, due to 



the short acceleration and cooling time-scales characterizing 
this system, especially close to the periastron ( |Khangulyan| 



et al. 20071. Observations show that this system is indeed 



a strong source of nonthermal time-dependent emission ex- 



tending from radio to TeV gamma-rays (see e.g. |Neronov fc| 



Chernyakova ( 2007a Q. 



Two variable TeV galactic gamma-ray sources, LS 5039 
and LSI 61 303 (see e.g. |Paredes| pOOo) ), are discussed as 
possible, although less evident, candidates representing "bi- 
nary pulsar" source population ( |Dubus|2006} |Mirabel|2006 



|Neronov fc Chernyakova|2007b| ). LS 5039 consists of 06.5V 
star and an unidentified compact object in a 3.9 day orbit. 
This object has been detected as the source of gamma rays 
by EGRET (so urce 3EG J1824-1514) JParedes et al. |2000[ ) 
and by HESS ( | Aharonian et al. ||2005| . LSI 61 303 is a bi- 
nary system with a BOVe star in a 26.5 day orbit. This source 
presumably associates with a low energy gamma-ray source 



2CG 135+01/3EG J0241+6103 (Maraschi fc Treves 1981 



Tavani et al. 19981. Recently, it has been detected in TeV 
gamma-rays ( Albert et al. |[2006[ ). The nature of the com- 
pact objects (black hole or a neutron star) in both sources 
is not yet firmly established. 

The collision of supersonic winds from two stars located 
in a binary system results in the formation of two terminat- 
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ing shock fronts and a tangential discontinuity separating 
relativistic and nonrelativistic parts of the shocked flow. It 
is believed that the shocked flow should propagate into a 
limited solid angle with rather high velocity. However the 
properties of the flow downstream the shock is unknown. 

The collision of the pulsar wind with the supersonic flow 
of the nonrelativistic plasma has been recently numerically 
modeled by several groups ( Bucciantini 2002 , S waluw|2003 



Vigelius et al. 2006 Romero et al. 20071 using different 



nonrelativistic versions of hydrodynamical codes. However, 
the dynamics of the relativistic plasma is rather specific, 
and the use of nonrelativistic codes cannot be a priori jus- 
tified. The modeling of collision of a relativistic pulsar wind 
with a nonrelativistic one demands an adequate treatment of 
distinct features of relativistic outflows. Moreover, in many 
previous studies the stellar wind has been approximated as 
plane parallel, while both winds initially expand radially. 
However, this purely geometric difference in the formulation 
of the problem, leads, in fact, to significantly different results 
and conclusions. 

In this paper we present the results of our studies con- 
ducted in the hydrodynamical limit, i.e. the role of the mag- 
netic field in dynamics has been ignored. The impact of the 
magnetic field will be published elswhere. 



2 PROPERTIES OF THE PULSAR AND 

STELLAR WINDS OF PSR 1259-63/SS2883 

Below we describe the properties of PSR 1259-63/SS2883 in 
which the interaction of the pulsar and stellar winds leads to 
the nonthermal emission observed in radio, X-ray and TeV 
gamma-ray energy bands. 



2.1 Pulsar wind 

The rotational losses of pulsars are released in the form of 
electromagnetic and kinetic energy fluxes of the wind. In- 
terpretation of observations usually results in a conclusion 
that the energy flux in the wind at large distances from the 
pulsar is concentrated in the particle kinetic energy. The ra- 
tio of the electromagnetic energy flux over the kinetic energy 
flux is described using cr-parameter. Typically a is much less 
than 1. Therefore, it is natural in the first approximation to 
consider the pulsar wind as purely hydrodynamical. Below 
we assume that the wind is ejected isotropically with the 
total energy flux equal to the total rotational losses of the 
pulsar E ro t- Thus the momentum flux density of the wind 
varies with the distance r to the pulsar as 



L = 



Erot 

4ncr 2 



(1) 



For a given Lorentz factor of the wind To, the proper density 
of the plasma in the pre-shock region is equal to 



En 



4irmc s rh 



(2) 



2.2 Parameters of the stellar wind 

In the system PSR 1259-63/SS2883 the 47.8 ms radio pulsar 
rotates around the Be star in an elliptic orbit with a period 



3.4 yr. The eccentricity of the orbit is 0.87. The distance be- 
tween stars in apastron is ~ 10 14 cm (Johnston et al, 1994). 
The Be star has a mass ~ lOM© and radius R* — 6Rq. 
Due to the very fast rotation of the star, the mass outflow 
is strongly anisotropic. The mass flux is concentrated along 
the equatorial plane forming a disk-like flow. In addition, 
there is also an isotropic component of the wind (so-called 
polar wind) with smaller mass flux, but higher velocity. 

The disks in Be stars are dense and expand rather 
slowly. Typical velocities of the plasma in the disk are 150 — 
300 km/s (Wat ers et al.|1988 l. The velocity of the plasma in 
the polar wind is much higher achieving 1500 — 2000 km/s 
( Snow |1993"] ). The ratio of the mass flux density in the disk 
over the polar wind is estimated ( Waters et al.|1988 1 as 



ridVd 



30 - 100 



(3) 



The mass flux of the disk-like wind from Be SS2883 has 
been estimated by several authors. The estimates were made 
basically using eclipses by the disk-like flow of the pulsed 



radio emission from the pulsar. Following Johnston ( 1999 1 
the mass loss rate M = 5 x 1O _8 M0 yr -1 and the opening 
angle of the disk 8 ~ 15° . Then the ratio of the momentum 
flux-densities from the Be star and the pulsar, 77, in the case 
of pulsar-disk interaction, is 



V = 



E rot sin(6)/2) 
cMv 



(4) 



For the given parameters characterizing PSR 1259- 
63/SS2883, E mt = 8 x lO^ergs" 1 , 6 = 15°, and the ve- 
locity of the equatorial outflow v — 200 km/s, one obtains 
77 ~ 5.5 x 10~ 2 with an uncertainty by a factor of 2-3. In 
the polar wind the velocity is much larger, v ~ 2000 km/s, 
but, at the same time the mass-flux is reduced by a factor 
of 30-100. Therefore, in the polar wind of the Be star, 77 ap- 
pears in the range 0.2 - 0.6. Thus, given the uncertainties 
in the mass-flux from the Be star and the variation of the 
momentum flux with distance due to the acceleration of the 
stellar wind by radiation pressure, the parameter r\ along 
the pulsar orbit may vary between 10 -2 and 1. 



3 COMPUTATIONAL METHODS 

The scheme of the plasma flow formed at the collision of 
the pulsar wind with the wind from the Be star is shown 
in Fig. [T] Both winds are terminated with formation of 
two shock waves. The post shock flow of the nonrelativis- 
tic plasma is separated from the shocked relativistic plasma 
flow by a contact discontinuity. Within the hydrodynamical 
classification, the contact discontinuity is treated as tangen- 
tial discontinuity ( Landau Sz Lifshits |1987 l. 

A remarkable feature of the post-shock flow arising from 
the collision of the two supersonic winds is that the flow 
is subsonic only in the region close to the line connecting 
the stars (i.e. the symmetry line). Propagating downstream 
the flow crosses sound line (shown in Fig. [T] by the thick 
black lines) and becomes supersonic. We note that the loca- 
tions of the sound lines in the nonrelativistic and relativistic 
flows are different. While the flow in the subsonic region is 
described by elliptic-type equations, in the supersonic re- 
gion the flow is hyperbolic. Thus the plasma flow in the 
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of the separatrix characteristics. This extention is based on 
the fact that in the hyperbolic region a Cauchy problem 
can be formulated with the initial data specified at some 
line located slightly downstream of the separatrix charac- 
teristic. Thus, the solution of the problem of the mixed-type 
flow can be reduced to a two-step treatment. At the first 
step (nearest zone solution) the solution can be found in the 
limited computational domain using the relaxation method. 
At the second step (far zone solution) the solution for large 
distances downstream of the separatrix characteristics can 
be obtained using the approach of Bogoval ov fc Tsinganos| 
(1999]). 



Figure 1. The schematic representation of the flow structure at 
the collision of the winds. The red lines show the termination 
shock waves. The thin blue line shows the contact discontinuity. 
The green lines show typical flow lines in the post shock region. 
The thick blue line - the starting surface for the extention of the 
solution in the hyperbolic region. The thick black lines - sound 
surfaces. The red lines attached to them - separatrix characteris- 
tics. 



post shock region appears to be of mixed-type. The criti- 
cal lines or separatrix characteristics play a crucial role in 



the mixed-type hydrodynamical (Landau fc Lifshits ||1987 



Hayes fc Probstein|1959| and m agnetohydrodynamicai~i |Bo-| 
govalov||1994| Tsinganos et al. 19961 flows. The separatrix 
characteristics are located in the hyperbolic region, and in 
Fig. [l] they are shown by thin red lines. The separatrix char- 
acteristics and sound lines may coincide only in a particular 
case when they are orthogonal to the flow lines. The separa- 
trix characteristics separate the flow regions with different 
causal connection. In the "subsonic" region limited by the 
separatrix characteristics every couple of points are causally 
connected through hydrodynamical or MHD signals. Down- 
stream the separatrix characteristics any two points are 
causally disconnected unless they are located in the cone 
formed by two different characteristics. Moreover, in a case 
of steady state mixed type flow, the solution setting on the 
separatrix characteristics specifies a unique solution in the 
hyperbolic region. 

The most general numerical method of mixed-type 
problem solution is the relaxation method. In this method, 
the time dependent problem of the plasma flow is solved 
starting from some initial state. The method was used for the 
mixed-type equation solution by many groups in the context 
of the problem of astrophysical jet formation ( |Bogovalov fc 
Tsinganosl [1999] [Rornanova et al. ||1993| |Ouyed fc Pudrits" 



1997|[Krasnopolsky, Li fc Blandford 2000 1. However, the ap- 



plicability of this method is rather limited. Usually transient 
equations are solved in a rather small computational do- 
main. Typically, the outer boundaries of the domain are lo- 
cated in the region of the hyperbolic flow. Theoretically they 
can be located at arbitrarily large distances downstream of 
the critical surfaces. However, in practice they are located 
not very far away due to the limitations imposed by mem- 
ory and power of computers. At the same time, one is often 
interested in the plasma flow rather far downstream from 
the critical lines. Therefore, it was proposed by Bogovalov 
|fc Tsinganos| ( |1999[ ) to extend the solution obtained by the 
relaxation method to arbitrarily large distances downstream 



3.1 Method for the solution in the nearest zone 

The computational method used in this work was devel- 
oped for the solution for the problem of the interaction of 
the relativistic pulsar wind with the interstellar medium 
(|Koldoba, Kuznetsov fc Usty ugova 2002 Bogovalov et al. 
|2005p . This method has certain advantages compared to 
other approaches published in the literature. Namely: 
(i) The equations which describe the dynamics of the plasma 
are solved only in the shocked region. There is no need to 
compute the plasma flow in the pre-shock region as it is done 



in previous studies (e.g. Bucciantini (20021; Vigelius et al 



|] ( [2006] ); |Swaluw| ( [2003] )) because the flow in the pre-shock 
region is a priory known; 

(ii) All discontinuities are considered as having zero thick- 
ness mathematical breaks. This allows us to avoid numer- 
ical diffusion of discontinuities. This property is especially 
important for calculations of radiation from the shocked ma- 
terial; 

(iii) The method has no limitations on the Lorentz factor 



of the wind (e.g. in the approach suggested by Komissarov 
|fc Lyubarsky| ( |2003[ ) the Lorenz factor of the pulsar wind 
cannot exceed 10). 

We assume that the winds from the pulsar and Be star 
are isotropic. In this case flow formed after the interaction 
is axially isotropic. In this paper we assume that both winds 
are cold. Strictly speaking, the stellar wind is hot. However, 
the sound speed in the stellar wind is significantly below the 
bulk motion speed, making this approximation well justified. 

The relativistic hydrodynamical (RHD) equations de- 
scribing the post-shock flow are 



dnj Id 9 
at r or oz 







dy wv r 1 
dt r dr 



9 I 2 2 \ & 2 P 
r Iwy v r + p) + -r^w-y v r v z = - 
v ' oz r 



dt 



19 2 ,9 I 2 2 , \ n 

-rwy v r v z + 7T— (ury v z +p) = 
oz y 



r or 



(5) 
(6) 
(7) 
(8) 



— h - 7T™7 V r + 7T W 7 V z = . 

at r or oz 

Here p, w,n and v are the pressure, enthalpy, particle density, 
and flow velocity, respectively At the termination shock 
front the bulk motion energy is transformed into the energy 



1 Here all thermodynamic quantities have their proper values, p 
and w are taken per unit volume in the local rest frame. 
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of the plasma's chaotic motion. The efficiency of this trans- 
formation depends on the incident angle with which plasma 
crosses the shock wave. 

The plasma flow in the nonrelativistic region is de- 
scribed by the following equations: 

dpv z 



dp 1 drpvr 
dt r dr 



+ 



dz 







dpv 



1 d 

H r 

dt r dr 



pv, 



d 



dpv z 1 d 
dt r dr 



+ 



(pvl + p) = 



(9) 
(10) 

(11) 



d_ 
dt 



pv 



1 d 
r dr 

d 

dz 



pv 



e+P + 



pv 



+e+p =0 



(12) 



Here p and e are the densities of mass and thermal energy, 
respectively, and p (= 2/3e) is the pressure. 

Generally, one faces certain difficulties while solving nu- 
merically the equations which describe interaction of rel- 
ativistic and nonrelativistic winds. The post shock flow is 
described by the systems of equations which are essentially 
different for the relativistic and nonrelativistic flow domains. 
The locations of discontinuities (shock waves and the contact 
discontinuity) are unknown a priori. In addition, there are 
other difficulties related to large differences in velocities of 
the relativistic and nonrelativistic flows. The velocity of rela- 
tivistic plasma is equal to c in the pre-shock region and larger 
than c/3 = 10 10 cm/s in the post shock region. The veloci- 
ties of the nonrelativistic stellar wind in PSR 1259/SS2883 
are ~ 3 • 10 8 cm/s and ~ 10 s cm/s in the pre-shock and 
post-shock regions, respectively. The spacial scales in the 
nonrelativistic and relativistic flows are similar and have to 
be meshed with cells of similar space resolution. This means 
that at integration of the full system of equations with ex- 
plicit numerical schemes, the time step will be limited by 
the Courante condition in the region of the relativistic flow. 
In the region of the nonrelativistic plasma flow this time- 
step will be less than the Courant limited time-step by two 
orders of magnitude. 

To overcome all these difficulties, a special numerical 
method has been developed for integration of transient HD 
and RHD equations. The integration of the HD equations 
was performed in the region limited by the nonrelativistic 
termination shock wave and by the contact discontinuity. 
The RHD equations have been integrated in the region lo- 
cated between the relativistic termination shock wave and 
the contact discontinuity. The plasma flow outside these re- 
gions is known. 

An adaptive computational mesh has been used. The 
mesh boundaries were located on the shock waves and the 
contact discontinuity. The boundary node location varied 
with change of the discontinuity locations. Usually the nodes 
were located on the system of beams which were fixed. The 
nodes moved only along these beams. The beams were di- 
rected along an axis of symmetry of the problem or spread 
radially from a fixed point (pulsar position for example) as 
it is shown in Fig. [2] The nodes on the beams located on 
the discontinuities were specified. The internal nodes were 
distributed between them uniformly. 




Figure 2. Typical mesh in the computational domain of the near- 
est zone region. Only every tenth cell is shown. 



The numerical integration of the HD and RHD equa- 
tions on the adaptive mesh has been performed with devel- 
oped numerical scheme of Godunov type of the second order 
for spacial variables ( Godunov |1993 1 . The fluxes of the con- 
servative variables on the cell faces were calculated with a 
help of the approximate solution of the Riemann problem 
with the limitation of antidiffusion terms. In the nonrela- 



tivistic domain the Roe linearization ( Roe 1986 \ has been 
used. 

The location of the shock fronts and the contact discon- 
tinuity have been defined in the process of the solution. The 
discontinuity decay has been explored for any discontinuity 
to find the location and velocity at the next time step. The 
Riemann problem was solved exactly at the nonrelativis- 
tic shock position. The Riemann problem at the relativis- 
tic shock front was solved in the acoustic approximation. 
Schematically the decay of the discontinuity between the 
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Figure 3. The scheme of the RHD discontinuity decay. The dis- 
continuity can decay into the shock terminating the infalling cold 
plasma, tangential discontinuity and a shock or rearefication wave 
propagating in the shocked material. The case when the shock 
propagates in the shocked material is shown. In the linearized 
case there is no difference between the shock and the rearefica- 
tion wave in the analytical representation. 



cold ( p = ) relativistic plasma inflow (state "0") and hot 
relativistic plasma outflow (state " 1" ) is shown in Fig. [3] At 
the calculation of the decay problem it was assumed that: 

(i) the termination shock is close to the steady state and the 
variation of all parameters can be described by linearized ( 
in relation to the steady state) Gugonio relationships. 

(ii) the variation of the parameters in the post shock region 
behind any rarefication or shock wave can be described by 
the linearized equations of RHD. 

Amplitudes of the waves propagating in the states " 0" 
and the state "1" are taken to satisfy the conditions at 
the contact (tangential) discontinuity: pressure and the nor- 
mal component of the velocity are constant. The diagram 
of the RHD discontinuity decay in variables (p, v) is shown 
in Fig. [4] The curve a is the Taube shock adiabata, go- 
ing from the point A - initial state of the plasma in the 
pre-shock region. Point E corresponds to the case when the 
termination shock wave is in the steady state. In the last 
case p = (2Q/3) cos 2 down stream the shock, where Q is 
energy density in the ultrarelativistic wind, ^ is the angle 
between the flow line and normal to the shock front. The 
curve j3 is the shock adiabata going from the point B - the 
state of the plasma downstream the discontinuity. The point 
C, where the curve a crosses the curve (5 corresponds to the 
state at the perfect solution of the Riemann problem. In our 
approach the state between the waves corresponds to the 
point D where the tangents to shock adiabates in the points 
E and B cross each other. It follows from the figure that the 
closer the points E and B placed to each other more accurate 
are the solution. 

Nonrelativistic HD equations allow a simple scaling of 
variables: t — > t/a, v — > av, p — > p/a 2 , p —* p. Moreover, 
this scaling does not change the conditions [p] = 0, v n — 
(which are viable in the steady state regime) at the contact 
discontinuity. Thus, the difficulties related to large difference 
in velocities of the relativistic and nonrelativistic flows may 
be overcome by rescaling of initial conditions for Eqs.Q- 
(121 in such a way that v — > c in the nonrelativistic wind. 
This approach allows us also to stabilize the flow at the con- 
tact discontinuity which in fact is a tangential discontinuity. 




xcosy 



cos\|/ 



Figure 4. Diagram of RHD discontinuity decay in variables 



3.2 Parameterization of the solution 

The steady state equation systems which describe both 
the nonrelativistic and relativistic outflows can be scaled 
through the value of the momentum flux of the nonrela- 
tivistic wind, pv 2 , at the distance D between the star and 
the pulsar. Also, it is natural to normalize all geometrical 
variables to D. In all figures below all the distances are given 
in this normalization unless stated otherwise. With such a 
normalization pv 2 — 1 at the dimensionless distance D = 1. 
At the same time the momentum flux of the relativistic wind 
at D = 1 becomes equal to the parameter r\ given by Eq. Q . 
In this way the interaction of the two flows is described by 
two parameters: r\ and the initial Lorentz factor of the pul- 
sar wind To. In the limit To 3> 1, the enthalpy per particle 
w/n ^> 1, therefore the properties of the flow only weakly 
depend on the initial Lorentz factor 70, thus in the case 
of pulsar winds with To 2> 1, the ratio of the momentum 
fluxes of the two winds rj becomes the key parameter which 
describes the system of interacting winds. Finally we note 
that the position (the distance to the pulsar) of the contact 
discontinuity on the symmetry axis may be approximately 
defined as 



(1 + v^) 



(13) 



3.3 Method for the solution in the far zone 

For magnetohydrodynamic nonrelativistic and relativistic 
flows, a method (hereafter BT method) has been suggested 
by |Bogovalov fc Tsinganos ( 1999} ). Below we discuss a pure 
hydrodynamical version of this method which can be applied 
to the far zone. 

The flow in the far zone is supersonic. The system of 
equations describing the steady state flow is hyperbolic. 
Therefore, an initial-value Cauchy problem can be formu- 
lated for this flow. 

It is convenient to introduce the flux function ip as fol- 
lows 



where <p is the azimuthal unit vector. 



(14) 
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The stationary HD equations accommodate the energy 
conservation law 



(e + 7r)7 = W(il>) , 



(15) 



where e and n are the thermal energy and pressure per par- 
ticle. According to this law, the average energy of particles 
remains constant along the flow line. This equation should 
be supplemented by the relativistic relationship between the 
components of the four-velocity u 



7 



(16) 



For analysing the behavior of the plasma at large dis- 
tances, it is convenient to deal with this equation in an or- 
thogonal curvilinear coordinate system (ip, a) formed by the 
flow lines and by the lines perpendicular to them. While ip 
varies with the motion across the flow lines, the coordinate a 
varies with the motion along the flow lines. The geometrical 
interval in these coordinates can be expressed as 



2 7/2 i 2 i 2 , 2 i 2 

g^dip + g a da + r dip , 



(17) 



where , g\ are the corresponding line elements (compo- 
nents of the metric tensor) . 



According to Landau & Lifshits ( 1994 1 the equation 



= (where T 1] is the energy-momentum tensor and 
;k" denotes covariant differentiation in these coordinates) 



can be written in the form 

dp 2 1 dg a 

— - u w — 

dtp y g a dip 



0. 



(18) 



The unknown variables here are z(a,ip) and r(a,ip). 
The metric coefficient g a can be obtained from Eq. (181, 



g a = exp 



where 



G(a,i,) = 



G(a,ip)dip 



1 dp 



tpW dip 



(19) 



(20) 



The lower limit of the integration in Eq.(19l is chosen to 



be such that the coordinate a is uniquely defined. In this 
way a coincides with the coordinate z where the surface of 
constant a crosses the axis of rotation. 



The metric coefficient g^, can be obtained from Eq.( 141 



in terms of the magnitude of the poloidal velocity as follows 
' (21) 

(22) 
(23) 
(24) 

dr/da, z a — dz/da, 



9i> = • 

rnup 

The orthogonality condition 

^aXip ~t~ Z^Z^p — , 

and the relationships 



2 _ 2 . 2 

g a v a -f- z a 
results in 

Ztpga 



</</> 



and 



and 



2 



2 . 2 



ri,g a 



9i> 



with g a given by Eq.([19f. Here r a 

r-ip — dr/dip, z^ — dz/dip. 



Eqs.(24l should be supplemented by appropriate 



'0 



;. ! o 
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Figure 5. The position of the starting line (green line), flow lines 
forming the coordinate system, contact discontinuity (blue line) 
and shock fronts at 77 = 2.5 X 10 — 3 . Only every 10th flow line 
is shown here. Pulsar is located at the point with coordinates 
r = 0, z = 0. Be star is located at the point with coordinates 
r = 0, z = 20. Geometry is normalized by the distance between 
the pulsar and the terminating relativistic shock wave. 



boundary conditions and initial values on some initial sur- 
face of constant a located in the nearest zone, but down- 
stream from all the critical surfaces. The form of the initial 
surface of constant a was obtained numerically via integra- 
tion of the following equations 

dr nu z dz 



r(nup) 



dip 



r(nu p ) 



(25) 



where nu z and nu r as well as the integral W(tp) are given 
by solution in the nearest zone. 

The boundary conditions on the axis of rotation and 
the equatorial plane are the same as the conditions in the 
nearest zone. No conditions at infinity are specified. 

Fig. [5] demonstrates a typical structure of the solution 
and the coordinate system obtained with the BT method. 
The terminating shocks are well reproduced. The nonrel- 
ativistic and relativistic flow are perfectly separated. The 
distribution of the selected flow lines and correspondingly 
of the coordinate system lines was chosen to reproduce ac- 
curately the structure of the shocks. 



4 RESULTS 

The modeling of hydrodynamical collision of two winds from 
a massive star and pulsar has been performed for a wide 
range of r) parameter: from r\ = 10 -4 , to 77 = 20, which 
covers both cases when momentum flux is dominated by the 
wind of pulsar (77 > 1) or by optical star (77 < 1). 

Two computational methods have been applied for so- 
lutions in the near and far zones. The comparison of two 
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Figure 6. Comparison of the results obtained by different meth- 
ods in the general computational domain for r] = 1.110 - 3 . The 
color represents the distribution of temperature in the post shock 
region. The left panel shows the solution obtained by the relax- 
ation method (nearest zone solution). The middle panel shows 
the solution obtained in the far zone and the right panel shows 
the comparison of the solutions obtained by these methods. The 
solid lines here show the position of the flow lines and shock fronts 
obtained by the relaxation method. Color represents the solution 
obtained by BT method. In the region of the relativistic flow the 
temperature is expressed in the units mc 2 UQ. In the nonrelativis- 
tic region the temperature is expressed in the units mtij, where 
vq is the initial velocity. 



methods in the far zone is shown in Fig. [5] Although the two 
solutions are in good agreement, being practically identical, 
nevertheless there is a deference between these solutions, 
namely the shock structures predicted by two methods are 
somewhat different at the backside point of the pulsar wind 
termination shock. In the backside point the post-shock su- 
personic relativistic flow converges with the axis, thus a re- 
flecting shock is formed. The structure of the reflecting shock 
depends on the angle to the symmetry axis. This behavior 
of the flow is analogous to the reflection of the shock wave 



from a solid wall (Landau & Lifshits 11987). Two different 



structures are possible at the reflection: (i) perfect reflection 
and (ii) Mach reflection. In the last case a triple configura- 
tion of the shock front with a tangential discontinuity is 
formed (see fig[7|. In this structure a region of subsonic flow 
is formed. While the BT method for the problem solution 
in the far zone is applicable only for pure supersonic flow, 
the subsonic region is not reproduced properly. However, the 
subsonic domain occupies a small region on plots shown in 
Fig. |6j otherwise the BT solution correctly reproduces all 
features relevant to the interpretation of observations. 

The interactions of the pulsar wind with a plane par- 
allel nonrelativistic flow always result in the formation of a 
closed termination shock wave as shown in Fig. [8] Similar 
structure has been revealed in calculations by Bucciantini| 



| ( |2002[ ); |Vigelius et al~| ( |2006| ); |Swaluw| ( |2003[ ). However, 
this assumption cannot be applied to the radially expand- 
ing winds. The shock wave associated with the termination 
of the pulsar wind becomes closed only in the case when 
the momentum flux in the nonrelativistic wind significantly 
exceeds the momentum flux in the pulsar wind, namely for 
r\ < 1.25 ■ 10 -2 . For larger r\ the shock is open as shown in 
Fig-i 

The properties of the post shock flow can be character- 
ized by the location of the termination shocks at the sym- 
metry line and by the asymptotical properties of the shocks. 





Figure 7. The possible structure of the shock waves at the re- 
flection from the symmetry axis. The thick solid line - the shock 
fronts. The thin solid line - flow line. Dashed line - the tangential 
discontinuity. Two variants are possible. On the left panel the 
perfect reflection. On the right panel - Mah reflection. 



Dependence of the location of the shocks and the contact 
discontinuity on r\ at the symmetry line (at the bow shock) 
is shown in Fig. 1 1 1 1 At small rj, the position of the discon- 
tinuity relative to the pulsar rd ~ ^/rj- The distance to the 
backward point of the relativistic shock front does not follow 
this law; it diverges at r\ ~ 1.25 • 10 -2 . 

An important parameter characterizing the post-shock 
flow is the asymptotic opening angle 9. Dependence of 8 on 
77 is shown in Fig. [12] for calculated positions of the non- 
relativistic bow shock, the tangential contact discontinuity, 
and the relativistic bow shock. The function 6(rf) can be 
interpolated as 

6 = exp(4.7516)?7 ' 217 = 115.7693?? ' 217 (26) 



for the nonrelativistic shock wave, 
6» = 41.068 lg 77 + 71.693 
for the relativistic shock wave, and 
= 28.64(2 - t? 275 )^ 173 



(27) 



(28) 



for the contact discontinuity. 

Interestingly, Eq. [28] agrees quite well with the analyt- 
ical interpolation of the asymptotic opening angle of the 



contact discontinuity proposed by Eichler & Usov (19931, 
despite the fact that the interpolation of |Eichler fc Usov 



(|1993 1 is based on the results of simulations of Girard & 



Willson (1987) performed for collisions of nonrelativistic 



winds. 

It follows from Fig. [9] that the post-shock flow is accel- 
erated to rather large Lorentz factors. Even in the nearest 
zone (limited by the radial distance equal to the distance 
between the two stars from the symmetry axis) the bulk 
motion Lorentz factor increases up to 7 = 2. In the far zone 
the Lorentz factor becomes very large, e.g. in the zone lim- 
ited by r = 50 (the size of the computational domain) the 
bulk Lorentz factor can be as large as 100. The results of cal- 
culations shown in Fig. [9] are performed for 77 = 1 which for 
the system PSR 1250/SS2983 approximately corresponds to 
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Figure 8. The structure of the post shock flow in the nearest 
zone solution at the interaction of the relativistic wind with the 
plane parallel flow of the nonrelativistic plasma. The color rep- 
resents the Lorentz factor. Red thick lines show the sound lines. 
Geometry is normalized on the distance between the pulsar and 
the relativistic shock front at the symmetry axis. 



the case when the pulsar interacts with the polar wind of 
the Be star. In Fig.[lO]we show also the results for r\ = 0.05 
which corresponds to the location of the pulsar in the equa- 
torial wind of the Be star. This case is shown in Fig. |10| 
Even in this case, the post-shock flow can be accelerated in 
the downstream region to very large Lorentz factors. 

This interesting effect obviously is related to adiabatic 



losses. According to (151 the full energy per particle, which 



is the product of the bulk motion Lorentz factor 7 and en- 
thalpy (e + 7r), is conserved along the flow line. Thus the ac- 
celeration of the flow is reduced to the transfer of the thermal 
energy to the bulk motion. Initially the energy of particles in 
the cold pulsar wind is domined by the kinetic energy of the 
motion. At the shock the total energy per particle remains 
unchanged, but a part of the energy is transformed to the 
thermal energy of particles. 

The relationship between bulk motion Lorentz factor 




Figure 9. The post shock flow for r\ 
the bulk Lorentz factor. 



n=o.o5 



1. The color represents 




Figure 10. The post shock flow for rj = 0.05. The color represents 
the bulk Lorentz factor. 



and post shock enthalpy depends on the incident angle be- 
tween the shock front and the flow line. At the symmetry 
line, practically all bulk motion energy is transformed to the 
internal energy (enthalpy) of particles. Downstream of the 
shock the pressure reduces because of the flow expansion, 
and thus the bulk motion Lorentz factor increases. Corre- 
spondingly the temperature falls down, in other words we 
deal with adiabatic losses. Formally, the bulk motion Lorentz 
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Figure 11. rj dependence of the location of the shock waves at 
the symmetry line. The solid line corrcspoends to the analitical 
approximation defined by Eq . ( | 1 3 [ l ■ 



Figure 13. The post shocked flow in the far zone for r) = 1.2 ■ 
10 -4 . The color represents the Lorentz factor. 



t - relativistic shock 
• - contact discontinuity 
a - nonrelativistic shock 




Figure 12. T) dependence of the asymptotical opening angle of 
the shocks and discontinuities in the post shock flow together with 
interpolations defined by Eqs.(|26|-||28|l. 



factor can achieve the initial value provided that the post- 
shock is not terminated by the interstellar medium. 

The bulk motion acceleration to very large Lorentz fac- 
tors is possible not only at the collision of the pulsar and 
stellar winds of comparable power as shown in Fig. [9] Even 
in the case when the momentum flux in the Be wind strongly 
exceeds the momentum flux in the pulsar wind, the-post 
shock flow can be reaccelerated to 7 3> 1 as it shown in fig. 

LU 



5 DISCUSSION: IMPLICATIONS AND 
LIMITATIONS 

In this paper, we have conducted a detailed numerical 
study of the interaction of the relativistic and nonrelativis- 
tic winds assuming isotropic, radially expanding winds and 
ignoring the role of the magnetic field. In fact, the cold 
ultrarelativistic pulsar winds are generally believed to be 



highly anisotropic (Bogovalov 1999 Lyubarsky 20021. The 



anisotropy of the energy flux in the wind can result in a 
nonaxisymmetric form of the termination shock wave at the 
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vicinity of the symmetry axis, as it follows from the studies 
of Vigelius et al. (20061. Whether the anisotropy has a no- 



ticeable impact on the result and conclusions of this paper 
is a subject of further investigations and will be discussed 
elsewhere. This concerns also the role of the magnetic field. 
Generally, as it follows from calculations of the interaction 



of the pulsar wind with the interstellar medium (Bogov 



alov et al.|2005 , Komissarov & Lyubarsky 2003 ; Bucciantini 



|2002 I, the role of the magnetic field is rather important in 
the post shock region almost independent of the strength of 
the field. This is explained by fast amplification of the mag- 
netic field in the post shock region due to deceleration of the 
flow (Khanguly an fc Bogovalov| 2004[). However, as shown in 
this paper, the flow is not decelerated. Just the opposite - 
it can be accelerated to large bulk motion Lorentz factor. 
Thus although one should expect an impact of the magnetic 
field on the post shocked flow, this affect most likely will not 
be so strong as in the case of the interaction of the pulsar 
wind with interstellar medium. Detailed MHD calculations 
are needed to to clarify the role of the magnetic field. 

The effect of reacceleration of the post shock flow to rel- 
ativistic bulk motion Lorentz factors has direct implication 
to the interpretation of observations of high energy 7- and X- 
rays from binary pulsar systems like PSR 1259-63/SS2883. 
This effect strongly modifies the relationship between the 
synchrotron X-ray and inverse Compton gamma-ray fluxes 
produced by the same population of relativistic electrons. 
It is well known that in the pulsar wind nebulae (pleri- 
ons) which are formed, by the interaction of pulsar winds 
with the interstellar medium, the ratio between the X-ray 
and VHE 7-ray fluxes are defined by the ratio of the en- 
ergy density of the magnetic field to the energy density of 
soft radiation field, provided that Compton scattering takes 



place in the Thompson regime ( Aharonian, Atoyan fe Ki- 
fune |[l997 l. In binary systems inverse Compton scattering 
proceeds in the Klein-Nishina regime which changes the re- 
lationship between the X-ray and VHE gamma-ray fluxes 
( Khangulyan fc Aharonian|2005 1 . Due to our calculations it 
becomes clear that a significant deviation from the standard 
relations should be expected also from hydrodynamics. 

Let us consider these processes in more deatils assuming 
that the magnetic field is present in the wind. Relativistic 
particles moving in the magnetic field usually produce syn- 
chrotron radiation. However, if the wind is cold (in this case 
the velocity of particles coincides with the bulk motion ve- 
locity) these particles do not produce synchrotron radiation. 
However, they can produce gamma-radiation with a specific 



sharp spectral feature through the IC scattering ( Bogovalov 
fc Aharonian|[2000l) . An example of such a system is a cold 



pulsar wind which does not produce synchrotron radiation 
in the pre-shock region. In such a system the "standart" 
relation between synchrotron and IC radiation components 
is violated. Remarkably, the possibility for the formation of 
relativistic flows in the post-shock region, as revealed in this 
paper, shows that in binary pulsar systems we should ex- 
pect a "non-standard" relation between synchrotron and IC 
radiation components in post-shock region as well. 

Moreover, due to the large bulk motion Lorentz factor 
we should expect strong modulation of the observed non- 
thermal radiation of electrons. Indeed in the case of binary- 
pulsar systems the direction of the post shock flow varies 
with the motion of the pulsar along the orbit around the 



star. This implies significant changes of the Doppler fac- 
tor, 5, given the large value of the Lorentz factor. Namely, 
S <C 1 for large viewing angles <j> (e.g. close to 90°) and 
8 ^> 1 for small viewing angles. Correspondingly this will 
have a strong impact on the lightcurve of nonthermal radia- 
tion of electrons, F 7 oc S" where typically n ^ 3. This effect, 
in particular, can naturally explain the interesting feature 
of the nonthermal emission of PSR 1259/SS2883, the both 
synchrotron and inverse Compton components of which dis- 
appear during the periastron passage of the pulsar (see e.g. 
Neronov & Chernyakova ( 2007a l). This interesting issue will 
be discussed elsewhere. 
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